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Abstract 

Two procedures to compute the output distribution (ps of certain stack filters S (so 
called erosion-dilation cascades) are given. One rests on the disjunctive normal form of S 
and also yields the rank selection probabilities. The other is based on inclusion-exclusion 
and e.g. yields cjis for some important LC/i [/-operators S. Properties of (jis can be used to 
characterize smoothing properties. 



1 Introduction 



The LULU operators are well known in the nonlinear multiresolution analysis of sequences. 
The notation for the basic operators L„ and U^, where n E N is a parameter related to the 
window size, has given rise to the name LULU for the theory of these operators and their 
compositions. Since the time they were introduced nearly thirty year ago, while also being used 
in practical problems, they slowly led to the development of a new framework for characterizing, 
evaluating, comparing and designing nonlinear smoothers. This framework is based on concepts 
like idempotency, co-idempotency, trend preservation, total variation preservation, consistent 
decomposition. 



As opposed to the deterministic nature of the above properties, the focus of this paper is on 
properties of the LULU operators in the setting of random sequences. More precisely, this setting 
can be described as follows: Suppose that X is a bi-infinite sequence of random variables Xi 
{i G Z) which are independent and with a common (cumulative) distribution function Fx{t) 
from L^{[0, 1], [0, 1]). Let be a smoother. Then we consider the following two questions: 

1. Find a map cps ■ [0,1] [0,1] such that the common output distribution Fsxit) of {SX)i 
{i € I) equals Fsx = <Ps ° Fx- The function (ps is also called distribution transfer. 

2. Characterize the smoothing effect an operator S has on a random sequence X in terms of 
the properties of the common distribution of {SX)i (i E Z). 
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With regard to the first question we present a new technique which one may call "expansion 
calculus" which uses a shorthand notation for the probability of composite events and a set of 
rules for manipulation. Using this technique we provide new elegant proofs of the earlier results 
in [1] for the distribution transfer of the operators LnUn and (dually) UnLn- The power of this 
approach is further demonstrated by deriving the distribution transfer maps for the alternating 
sequential filters Cn = LnUnLn~iUn-i---LiUi and F„ = UnLnUn~iLn~i---UiLi. 

With regard to the second question, we may note that it is reasonable to expect that a smoother 
should reduce the standard deviation of a random sequence. Indeed, for simple distributions 
(e.g. uniform) and filters with small window size (three point average, Mi, LiUi, UiLi) when 
the computations can be carried out a significant reduction of the standard deviation is observed 
(for the uniform distribution the mentioned filters reduce the standard deviation respectively 
by factors of 3, 5/3, 1.293, 1.293). However, in general, obtaining such results is to a large 
extent practically impossible due to the technical complexity particularly when nonlinear filters 
are concerned. In this paper we propose a new concept which characterizes the probability of 
the occurrence of outliers rather then considering the standard deviation. We call it robustness. 
Upper robustness characterizes the probability of positive outliers while the lower robustness 
characterizes the probability of negative outliers. In general, the higher the order of robustness 
of a smoother the lower the probability of occurrence of outliers in the output sequence. In terms 
of this concept it is easy to characterize a smoother given its distribution transfer function. 

The paper is structured as follows. In the next section we give the definitions of the LULU 
operators with some fundamental properties. The concept of robustness is defined and studied 
Section 3. In Section 4 we demonstrate the application of inclusion-exclusion principle for deriv- 
ing the distribution transfer function of erosion-dilation cascades, a kind of operator frequently 
used in Mathematical Morphology. This method is considerably refined in Section 5 where it is 
applied to LULU-operators which are in fact particular cases of such cascades. Formulas for the 
major LULU-operators are obtained explicitly or recursively. Using these results the robustness 
of these operators is also analyzed. Section 6 proposes to substitute inclusion-exclusion by some 
principle of exclusion which is better suited for erosion-dilation cascades that don't allow the 
refinements of inclusion-exclusion possible for LC/LC/-operators. 



2 The basics of the LULU theory 

Given a bi-infinite sequence x = {xi)i^i and n G N the basic LULU operators L„ and Un are 
defined as follows 

{LnX)i := (Xi-n A Xi-n+1 A • • • A Xj) V (Xi-n+l A • • • A Xj+i) V • • • V (Xj A • • • A Xi+n) (1) 
{UnX)i := {Xi-n V Xi-n+1 V • • • V Xj) A {Xi-n+1 V • • • V Xj+i) A • • • A (Xj V • • • V Xi+n) (2) 

where a A (3 := min(a,/3), and a V /3 := max(a,/3) for all a,/3 S M. Central to the theory is 
the concept of separator, which we define below. For every a € Z the operator Ea : — )• 
given by 

{EaX)i = Xj+a, i € Z, 

is called a shift operator. 

Definition 1 An operator S : — > is called a separator if 
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(t) SoEa = EaoS, o G Z ; 



(horizontal shift invariance) 



(ti) P{f + c) = P{f) + c, f,c 



c- constant function (vertical shift invariance); 



(Hi) P{af) = aP{f), a G M, a > 0, / G M^; 



(scale invariance) 



(iv) PoP = P; 



(Idempotence) 



(v) {id -P)o {id -P) = id-P. 



( Co-idempotence ) 



The first two axioms in Definition [T] and partially the third one were first introduced as required 
properties of nonlinear smoothers by Mallows, [1]. Rohwer further made the concept of a 
smoother more precise by using the properties (i)-(iii) as a definition of this concept. The axiom 
(iv) is an essential requirement for what is called a morphological filter, [11], [12], [13]. In fact, a 
morphological filter is exactly an increasing operator which satisfies (iv). The co-idempotence 
axiom (v) in Definition [1] was introduced by Rohwer in where it is also shown that it is an 
essential requirement for operators extracting signal from a sequence. More precisely, axioms 
(iv) and (v) provide for consistent separation of noise from signal in the following sense: Having 
extracted a signal Sx from a sequence x, the additive residual {I — S)x, the noise, should contain 
no signal left, that is S o {I — S) = 0. Similarly, the signal Sx should contain no noise, that 
is {I — S) o S = 0. It was shown in [8j that L„, [/„ and their compositions LnUn,UnLn are 
separators. 

The smoothing effect of Ln on an input sequence is the removal of picks, while the smoothing 
effect of Un is the removal of pits. The composite effect of the two LU-operators LrJJn and UnLn 
is that the output sequence contains neither picks nor pits which will fit in the window of the 
operators. These are the so called n-monotone sequences, [8]. Let us recall that a sequence x is 
n-monotone if any subsequence of n+1 consecutive elements is monotone. For various technical 
reasons the analysis is typically restricted to the set Mi oi absolutely summable sequences. Let 
M.n denote the set of all sequences x G A^i which are n-monotone. Then 



is the set of signals. 

The power of the LC/-operators as separators is further demonstrated by their trend preservation 
properties. Let us recall, see [8], that an operator is called neighbor trend preserving if 
{Sx)i < {Sx)i+i whenever Xj < Xj+i, i G N. An operator S is fully trend preserving if both 
S and I — S are neighbor trend preserving. The operators L„, Un and all their compositions are 
fully trend preserving. With the total variation of a sequence. 



a generally accepted measure for the amount of contrast present, since it is a semi-norm on 
A^i, any separation may only increase the total variation. More precisely, for any operator 
S : Ml — )■ Ml we have 



Mn = Range{LnUn) = Range{UnLn) 
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(3) 



All operators S that are fully trend preserving have variation preservation, in that 





(4) 
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We mention these properties because they provide but few of the motivation for studying the 
robustness of operators, when the popular medians are optimal in that respect. We intend to 
show that some Lt/LCZ-composition are nearly as good as the medians, but have superiority 
most important aspects. 

An operator S satisfying property dU is called total variation preserving, [6]. As mentioned 
already, the LC/-operators are total variation preserving. 



3 Distribution transfer and degree of robustness of a smoother 

Suppose that X is a bi-infinite sequence of random variables Xi (i € Z) which are independent 
and with a common (cumulative) distribution function Fx- Let 5" be a smoother. As stated in 
the introduction we seek a function (ps '■ [0,1] — )• [0,1], called a distribution transfer function 
such that 

Fsx = 4>s°Fx (5) 

is the common distribution of {SX)i (i G Z). We should note that for an arbitrary smoother 
the existence of such a distribution transfer function is not obvious. However, for the smoothers 
typically considered in nonlinear signal processing (i.e. stack filters of which the LULU operators 
are particular cases) such a function does not only exist but it is a polynomial. For example, it 
is shown in [5j that the distribution transfer function of the ranked order operators 

{Rnkx)i = the kih. smallest value of 

is given by 

</'«..(p)= E V(i-p)'"+'-^'- (6) 

j=k V / 

The popular median smoothers M„, n G N, are particular cases of the ranked order operators, 
namely M„ = Rn,n+i- Hence we have 

j=n+l V / 

Note that in terms of ([5]) the common distribution function of (M„X)j, z G Z, is 



j=n+l 

Using that 



FMMt)= E )Fi,m-Fx{t)r^'-\t^^. 



^0Mjp) = (2n + l)n)p"(l-p)' 



its density is 



fM^t) = -ct>MSFx{t))fx{t) = {2n + l){ )FjJ(t)(l - Fx{t)rfx{t) 



4 



where fx{t) = ^Fx{t), t G M, is the common density of Xj, i S Z. The distribution of the 
output sequence of the basic smoothers L„ and C/„ is derived in [8|. Equivalently these results 
can be formulated in terms of distribution transfer. More precisely we have 



A primary aim of the processing of signals through nonlinear smoothers is the removal of im- 
pulsive noise. Therefore, the power of such a smoother can be characterized by how well it 
eliminates outliers in a random sequence. The concepts of robustness of a smoother introduced 
below are aimed at such characterization. 

Definition 2 A smoother S : — > is called lower robust of order r if there exists a 
constant a > such that for every bi-infinite sequence X of identically distributed random 
variables Xi (i £ there exists to £ such that P{Xi < t) < e implies P{{SX)i < t) < ae^ 
for all t < to and e > 0. 

Similarly, a smoother S : — )■ is called upper robust of order r if there exists a constant 
a > such that for every bi-infinity sequence X of identically distributed random variables Xi 
(i e 'Z) there exists to such that P{Xi > t) < e implies P{{SX)i > t) < ae^ for all t > to 
and e > 0. 

A smoother which is both lower robust of order r and upper robust of order r is called robust 
of order k. 

The reasoning behind these concepts is simple: If a distribution density is heavy tailed, there is 
a probability e that the size of a random variable is excessively large (larger than t) in absolute 
value. Using a non-linear smoother we would aim to restrict this to an acceptable probability 
ae^ that such an excessive value can appear in SX, by choosing a smoother with the order of 
robustness k. 

Clearly there is a general problem of smoothing: a trade-off to be made between making a 
smoother more robust, and the (inevitable) damage to the underlying signal preservation. (A 
smoother clearly cannot create information, but only selectively discard it.) This is fundamental. 
There are two main reasons for using one-sided robustness: Firstly, the unreasonable pulses often 
are only in one direction, as in the case of "glint" in signals reflected from objects with pieces 
of perfect reflectors, and there clearly are no reflections of negative intensity possible. Secondly, 
we may chose smoothers that are not symmetric, as are the LCZ-operators, for reasons that are 
of primary importance. In this case the robustness is determined from the sign of the impulse. 

The robustness of a smoother can be characterized through its distribution transfer function as 
stated in the theorem below. 

Theorem 3 Let the smoother S have a distribution transfer function (ps- Then 
a) S is lower robust of order r if and only if (psip) = 0{p'') as p — > 0. 




(9) 
(10) 
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b) S is upper robust of order r if and only if — p) — 1 = 0{p^) as p — > 0. 

Proof. Points a) and b) are proved using similar arguments. Hence we prove only a). Let 
4's{p) = 0{p^) as p — ?• 0. This means that there exists a > and 6 > such that (j)s{p) < ctp^ 
for all p G [0,(5). Let X be a sequence of identically distributed random variables with common 
distribution function Fx- Since limt^^oo Fx{t) = 0, there exists Iq such that Fx{to) < 6. Let 
t < to and e > be such that P{Xi < t) < e. The monotonicity of Fx imphes that Fx{t) G [0, 6). 
Then 

P{{SX)i <t) = Fsx{t) = MFx{t)) < a{Fx{t)f < as'', 

which proves that 5 is lower robust of order k. It is easy to see that the argument can be 
reversed so that the stated condition is also necessary. ■ 

In the common case when the distribution transfer function is a polynomial, conditions a) and 
b) can be formulated in a much simpler way as given in the next corollary. 

Corollary 4 Let the distribution transfer function of a smoother S be a polynomial (ps ■ Then 

a) S is lower robust of order r if and only if p = is a root of order r of <j)s- 

b) S is upper robust of order r if and only if p = 1 is a root of order r of cjjs — I- 

Using the distribution transfer functions given in ([9]) and (|10p it follows from Corollary [J] that 
Un is lower robust of order n + 1 and that L„ is upper robust of order n + 1. 

The robustness of the median filter Mn can be obtained from ([7]). Obviously p = is a root of 
order n+ 1. Furthermore, = 1. Then using also that p = 1 is a root of order n of -^(pMn: 

see ([8]), we obtain that p = 1 is a root of order n + 1 of (pM„ — 1- Therefore, Mn is robust of 
order n + 1. 

Clearly with symmetric smoothers, in that S{—x) = —S{x), the concepts of lower and upper 
robustness are not needed, as is the case for example with M„. However, we have to recall in 
this regard that the operators L„, C/„ and their compositions, which are the primary subject of 
our investigation, are not symmetric. A useful feature of the lower and upper robustness is that 
it can be induced through the point-wise defined partial order between the operators. Let us 
recall that given the maps A,B : — >■ M^, the relation A < B means that Ax < Bx for all 
X G R^. 

Theorem 5 Let ^,_B : be smoothers. If A < B then cpB < (pA- 

Proof. Let X be a sequence of independent random variables Xi {i G Z) uniformly distributed 
on [0, 1] . Let p G [0, 1]. If t is such that p = Fx{t) then 

Mp) = MFx{t)) = FBx{t) = P{{BX)i<t) 

< P{{AX),<t) = FAxit) = MFx{t)) = ct>A{p). 
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As a direct consequence of Theorem and Theorem [3] we obtain the following theorem. 

Theorem 6 Let A,B -.R^ be smoothers such that A< B. Then 

a) If A is lower robust to the order k, then so is B. 

b) If B is upper robust to the order k, then so is A. 

Using Theorem[6]one can derive statements about the lower robustness and the upper robustness 
of the LC/-operators: 

UnLn <Mn< LJJn (11) 

Therefore, UnLn inherits the upper-robustness of M„, while Lnlln inherits the lower-robustness 
of Mn- More precisely 

• UnLn is upper robust of order n + 1; (12) 

• LnUn is lower robust of order n + 1 

One may expect that, since L„ is upper robust of order n + 1 and Un is lower robust of order 
n + 1, their compositions should be both lower and upper robust of order n + 1. However, as 
we will see later, this is not the case. The problem is the following. The definition of robustness 
requires that the random variables in the sequence X are identically distributed but they are not 
necessarily independent. However, the distribution transfer functions (pi^ and <j)u„ are derived 
under the assumption of such independence. Noting that entries in the sequences LnX are not 
independent, it becomes clear that the common distribution of UnLnX cannot be obtained by 
applying to F^^x- More generally, since the distribution transfer functions are derived for 
sequences of independent identically distributed random variables the equality 4)ab = 4>a° 4>b 
does not hold for arbitrary operators A and B. Therefore the order of robustness of B is not 
necessarily preserved by the composition AB. 

Observe that another concept of robustness is introduced in [10]. Other than Definition 2 it 
only applies to stack filters. The concept is similar in that it also based on certain probabilities 
(in this case "selection probabilities"). 



4 The output distribution of arbitrary erosion-dilation cascades 

Here we present a method for obtaining output distributions of so called erosion-dilation cas- 
cades (defined below). It essentially uses the inclusion-exclusion principle for the probability of 
simultaneous events. For convenience we recall this principle below. For n = 2 the easy proof is 
given later on along the way. 

Lemma 7 For any random variables Zi, Z2 ■ ■ ■ Zn it holds that 
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n 

P{Zi,---,Zn<t) = l-Y^P(Z,>t)+ Yl P{Zi,Zj>t) 

- P{z,,Zj,Zk>t) + ■■■ + {-iYP{Zi,---,Zn>t) 

l<i<j <k<n 

Let us recall that in the general setting of Mathematical Morphology [12] the basic operators 
Ln and Un are morphological opening and closing respectively. As such they are compositions 
of an erosion and a dilation. More precisely, for a sequence x = {xi)i^i we have 

{Lnx)i := (V"(A"^)). 

iUnX)i := (A"(V"^))i 
where (A" x)^ := Xi-n A Xi-n+i A ■ ■ ■ A Xi 

is an erosion with structural element W = {—n, —n — 1, 1,0} and 

(V"^)j := XiV Xi+iV ■ ■ ■ W Xi+n 

is a dilation with structural element W' = {0, l,...,n}. Generalizing the -Lf7-operators LnUn 
and UfiLn, call a LULU- operator any composition of the basic smoothers L„ and Un, such as 
L^U/iL2UiU^. In particular, each L?7LC/-operator is a composition of dilations and erosions, 
that is, an cascade of dilations and erosions (CDE). More generally, each alternating sequential 
filter (ASF), which by definition is a composition of morphological openings and closings with 
structural elements of increasing size, is a CDE with the extra property of featuring the same 
number of erosions and dilations. See also [3j. 

We will demonstrate our method on two examples of cascades - the first in one dimension, the 
second in two dimensions. This method is considerably refined in the next section. 

Example 1. Consider S := V A V • I* 

cascade of an erosion f\ and dilations \l ,\l ■ 
It is not an ASF. To compute the distribution transfer of 5, let X be a bi-infinite sequence of 
independent identically distributed random variables X^. Put 

Y,= {yX^, Zr.= (^l\Y^, Ar.= {\IZ^. (13) 

Thus y, Z, A are again bi-infinite sequences of identically distributed (though dependent) random 
variables. Let t € M and p = Fx{t). Then 

4>s{p) = Fsx{t) = FA{t) = P{Ao<t) 

= P{Zo y Zi<t) = P{Zo < t and Zi < t) 

In order to reduce the Zj's to the l^'s we switch all < t to > t by using exclusion- inclusion (the 
case n = 2 in Lemma 7): 

P{Zo, Zi<t) = P{Zo <t)- P{Zo <t,Zi> t) 
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= P{Zo <t)-i P{Zi >t)- P{Zi, Zo>t)) 
= l-P{Zo>t)-P{Zi>t) + P{Zi,Zo>t) 

Since our Zj's are identically distributed we have P{Zo > t) = P{Zi > t) and hence 
(f>sip) = l-2P{Zo>t)+PiZi,Zo>t) 

= 1 - 2P(y_2 A y_i AYo>t) + p(y_i aYqa Yi, y_2 a y.i a > 
= 1 - 2P(y_2, Y_i,Yo >t) + p(Y_2, y_i, yo, n > t) 

= l-2PiYo,Y^,Y2>t) + P{Yo,Yi,Y2,Ys>t) 
By the dual of Lemma 7 and because e.g. P{Yo, Yi < t) = P{Yi,Y2 < t) = P{Y2, Y3 < t) we get 

(l^sip) = l-2{l-W{Yo<t) + 2P{Yo,Yi<t) + P{Yo,Y2<t)-P{Yo,Yi,Y2<t)) 

+ ( 1 - AP{Yq <t) + SP{Yo, Yi<t) + 2P{Yo, Y2<t) + P{Yo, Ys<t)- 2P{Yo, Yi,Y2 < t) 
-2P{Yo, Yi,Ys <t) + P{Yo, Yi,Y2, Ys<t)) 

= 2P{Yo <t)- P{Yo, Yi<t)+ P{Yo, Y3<t)- 2P{Yo, Yi,Ys <t)+ P{Yo, Yi,Y2,Y3 < t) 

= 2P{Xo, Xi,X2, Xs<t)- P{Xo, Xi, X2, X3, X4 < i) + ^'(^o, Xi, • • • , Xe < t) 
-2P(Xo, Xi, • • • < t) + P{Xo, Xi, • • • , Xe < 

= - + - 2p'^ + 

= 2p4-/ 



Example 2. Let S be an opening on M^^^ with defining structural element a 2 x 2 square. 

Let now X be an infinite 2-dimensional array of independent identically distributed random 
variables X^jj) where ranges over Z x Z. In order to derive the output distribution of S 
we put 
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Let t € M and p = Fx{t). The output distribution of 5 is 

= p(y(o,o) , 1^(1,0) , %-i) , < t) 

Following [1], which introduced that handy notation in the 1-dimensional case, we abbreviate 
the latter as 

((0,0), (1,0), (0,-1), (l,-l))y 
If say (((0,0),(l,0),(0,-l))y means 

P(>^(i,o) <i,%o),%-i) 

then it follows from Lemma 7 and from translation invariance (e.g. ((0, 0), (1, 0)) = ((0, —1), (1, — 1)) 
that 

Mp) = ((0,0),(l,0),(0,-l),(l,-l))y 

= 1 - 4(o:o)y + 2((o;o) , (T;o))y + 2((o;o) , jo~^)y + ((lyoj , jo~^)y 
+((o;o), (T~^)y - ((o;o), (o;^, (w))y - (m, (t;o), (i~^)y 
-(M, (o,-i),(T;o))y - ((i;o), (o,-i),(i,-i))y + ((o;o),(t;o), (o,-i),(i,-i))y 

According to the definition of ^(ij) we e.g. have 

((MJ, (0,-l))y = ( (0:0),(-l,0),(M),(-l,l), (0,-l),(-l,-l),(o;oy,(-l,0) )x 
= ((0;0), (-1,0), (oj), (-1, 1), (0, -1), (-1, -l))x 
Putting q = 1 — p = P(X(o,o) > t) the latter contributes a term to 

,/.l(p) = l-Aq^ + 2q^ + 2q^ ^q"^ + q'^ -q^ -q^ -q^ -q^ + 
= l-4:q'^ + Aq^ + 2q^ - + q^ 
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5 Formulas for the distribution transfer of the major LULU- 
operators 



As it was already done in the preceding section it is often convenient to use the notation q = 1—p. 
For example the output distribution of M„, L„ and C/„ given in ([7]), (l9|) and ([10]) respectively 
can be written in the following shorter form: 

<PmJp) 
<Pu„ip) 

Theorem 11 below deals with the output distribution of LnUn and UnLn- They were first derived 
in [2], but the statement of the theorem was also independently proved by Butler [1]. In 5.1 
we present a proof using Butler's "expansion calculus". In 5.2 this method is applied to more 
complicated situations. 

5.1 The output distribution of the Lf/-operators 

First, observe that instead of winding up with full blown inclusion-exclusion when switching 
all inequalities > t to < t (dual of Lemma 7), one can be economic and only switch some 
inequalities: 



(0,l,---,n)x = (0,---n-l)x-(0,---,n-l,n)x (17) 
= (0, • • • , n - 2)x - (0, • • • , n - 2, n - l)x - (0, • • ■ , n - 1, n)x 

= (0)x - (0, l)x - (0, T, 2)x (0, • • • , ^rn:, n)x 

n-l 

= l-(0)x-^(0,---,i,i + l)x 

i=0 



Lemma 8 [1, Corollary 4]: Let X be a bi-infinite sequence of random variables. Then 

n — 1 

{0,T,---,n)x = 1- [n + l](0)x + J^[n-i](0,T,---,i,i + l)x 

i=0 

Note that for i = we get the summand n(0, l)x- 



E (14) 

j=n+l V J / 

l-(n+l)g"+i+ng"+2, (15) 
(n+l)p"+i (16) 
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Proof: Prom 

{k + \)x = {k + l,k)x + {k + \,k)x 



= {k + l,k)x + {k + l,k,k-l)x + {k + l,k,k-l)x = 



= {k + l,k)x + {k + l,k,k-l)x + --- + {k + l,k,---,\,Q)x + {k + \,k,---, 0)x 
follows, by translation invar iance, that 

k-l 

i^,--- ,k,k + l)x = (k + l)x -{k,k + l)x - ^{k - i - l,k - i,- ■ ■ ,k,k + l)x 



i=0 



k-l 



= (0)x - (0, l)x - 5^(0, 1, • • • , z + 1, z + 2)x 

i=0 

Using (17) one derives for (say) n = 4 that 

3 

(0,T,2,3,4)x = l-{0)x-J2iO,---,k,k + l)x 



k=0 

3 



l-{0)x-Yl 



k-l 



{0)x - (0, l)x - ^(0, 1, • • • , i + 1, i + 2)x 



fe=0 

= l-5(0)x+4(0,l)x 



i=0 



+(0,1, 2)x 



+(0,l,2)x + (0,l,2,3)x 



+(0, 1, 2)x + (0, 1,2, 3)x + (0, 1,2,3, 4)x 



l-5(0)x + ^[4-i](0,l,---,i,z + l)x 

i=0 



□ 



Unsurprisingly, for dependently distributed random variables Bi certain combinations of Bj's 
being < t and simultaneously other Bj^s being > t, are impossible, i.e. have probability 0. More 
specifically: 

Lemma 9 [1, Theorem 10]: Let A he a bi-infinite identically distributed sequence of random 
variables and let B = \J^ A. Then 



(0, l,---,n- l,n) 



, n < r + 1 



(0, • • • ,r, r + 1, n — 1, n, ■ • • , n + r)A , r + l<n<2r + 4 
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For instance, for n = 5, r = 1 we have r + l<n<2r + 4, and so 

(0,T,2,3,4,5)b = (0,1,2,4,5,6U 

Let us give an ad hoc argument which conveys the spirit of the proof. In view of Bi = AiV Aj+i 
one e.g. has that < t <^ A^, Aq < t. Using inclusion-exclusion we get 

(0, 5, 2, 4, T, 3)b = (0, 5,2, 4)b - (0, 5, 2, 4, 1)b - (0, 5, 2, 4, 3)^ + (0, 5, 2, 4, 1, 3)^ 

= P{Ao, Ai,A5, Aq < t, B2,B4 >t)- P{Ao, Ai,A2, A5, Aq < t, B2, B4 > t) 

-P{AQ,Ai,As,A4,A5,AQ<t,B2,B4>t) + P{Ao,Ai,---Ae<t,B2,B4>t) 

Since A4, ^5 < t is incompatible with B4 = A4^V A^ > t, the last two terms are 0. Furthermore, 
given that A^ < t, the statement B^ > t amounts to A4 > t. Ditto, given that A2 < t, the 
statement B2 > t amounts to A3 > t. Hence 

(0,5,2,4,T,3)s 

= P{Ao,Ai,A5,Ae <t,A4> t,B2 > t) - P{Ao,Ai,A2,A5,Aq <t,A4> t,^3 > t) 
= ( P{Ao,Ai,A5,Aq <t,A4> t) - P{Ao,Ai,A5,Ae <t,A4> t,A2 < t,A3 <t) ) 

-P(^0, ^1,^,^6 <t,A4> t,A2 < t,A3>t) 

= P{Ao, Ai,A5, Ae <t,A4>t)- P{Ao, ^1,^5, Aq < i, ^4 > t, A2 < t) 

= (0,1,2,4,5,6)a 
Dualizing Lemma 9 yields: 

Lemma 10 [1, Corollary 11]: Let B = /\'' A. Then 

'0 , n<r+l 

< 

(0, • • • , r, r + l,n — l,n, • • • , n + r)^ , r + l<n<2r + 4 



(0, 1 ■ ■ ■ ,n - l,n)B = < 



Theorem 11 The distribution transfer functions of LnUn and UnLn are: 

H^uM = + np^+\ + p^^+\ + l(n - l)(n + 2)p'^+\\ (18) 
4>U^lAp) = 1 - ^L^uAq) = 1 - - npq^+' - pq^^+' - ^{n - l)(n + 2)p\^^+\ (19) 
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Proof: Since L„i7„ = (V" A")(A" V") = V" A'" V" we put 

n 2n n 

A = \/X, B = /\A, C = \Jb 

and calculate 

HnuM = P{Co<t) = (0)c = (0,---,n)B 

n-l 



= 1- [n + l](0)B + ^[n-i](0,l,---,z,i + l)B (dual of Lemma 8) 

n-l 

= 1 - [n + 1](0, • • -,2^) A + ^^(0, • • • ,2n + 1)a + XI ^ (Lemma 10, r = 2n) 



i=l 



= l-[n + l] 



2n-l 



1 - [2n + 1](0)a + X [2n - i](0, 1, • • • + 1)a 

j=0 

2n 

1 - [2n + 2] (0)a + X [2n + I - ^] (0, T, • • • J, i + 1)a 

i=0 

2n 

= [n + 1](0)^ + - ^](0 J, • • • J,i + 1)^ 

j=0 

n 

+ 1](0U - n(0, 1)a + ^[i - n](0,T, • • • + 1)^ 
j=i 

2n 

+(0,T,---,^Hn:,n + 2)A+ XI [«-?^](0,T,--- + 



(Lemma 8) 



= n 



i=n+2 



[n + 1] (0, • • • , n)x - n(0, • • • , ra + l)x + + (0, • • • , n, n + 1, n + 2, • • • , 2n + 2)x 



2?i 



+ X — ra](0, • • • ,n, n + l,i,i + 1, • • • ,i + 1 + n)x (Lemma 9) 

i=n+2 

2n 

= (n + l)p-+i-np"+2+p2-+2^+ X(^-^)^^'"^V 

i=ra+2 

= + np"+i - np"+^p + p2"+2g + (2 + 3 + . . . + n^^+^g^ 



□ 



From (18) it is clear that is the highest power of p dividing (I>l„u„- An easy calculation 

confirms that, as a polynomial in p, the right hand side of (19) is (2n + 3)p^ + (• • ■)p'^ + • • •• 
Prom Corollary 4 hence follows that LnUn is lower robust of order n + 1, but upper robust only 
of order 2. 
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5.2 The output distributions of the L?7LC/-operators C„ and F„ 

We consider next the specific, mutually dual LULU-operators 

Cn = LnUnLn-lUn-1 • • • LiUi, 



In view of 



J^n = UnLjiUn-lLn-l ■ ■ -UiLi. 



A V ^r^-2 

Vn . 2n. .n 
A V 

Vn A 2nv /2n— 1 » 2n— 2i /n— 1 
A V A V 



we define the following doubly infinite sequences of identically distributed random variables. 
Starting with a sequence X of i.i.d. random variables, put 



Vn— 1 

• 2n-2 
B := l\ A 

. /2n-l 

C := V ^ 

^ - A ^ 

Vn 



Theorem 12 /i, Theorem I4] With A,B as defined above the output distribution (l>Cn of Cn 
can be computed recursively as follows: 

= 0C„_i + n{G2n - G2„-l), 

where 

G2n ■■= (0,---,2n-l,2^,2n + l,-- -,471)5 



G2„-i := (0,---,2n-2,2n- l,2n,---,4n-2)^ 
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Proof: First, one calculates 

= (O)c' = (0,---,n-l)B 



n-2 



= 1 — n(0)s + — 1 — ■i](0, 1, • • • ,i, i + l)s (dual of Lemma 8) 
j=o 

n-2 

= 1 -n(0,---,2n-2)A + [n- l](0,---,2n- 1)a + 5^0 (Lemma 10, r = 2n - 2) 



1=1 



= 1 - n(0, ■ ■ ■ , 2n - 2)a + [n - 1](0, • • • , 2n - 1)a 
The expansion of (pCn is driven a bit further: 
<t>C^ = (0)E = (0,---,n)i, 



(20) 



This yields 

(l)C„ - nG2n 



n—1 



= 1 - [n + 1](0)£) + ^[n - i](0, 1,--- + (dual of Lemma 8) 

i=0 

n-l 

= 1 - [n + 1](0, • • • , 2n)c + n(0, • • • , 2n + l)c + ^ (Lemma 10, r = 2n) 



i=l 



1 - [n + 1] 



2n-l 



+n 



1- [2n + l](0)c+ ^[2n-z](0,l,---,i,z + l)c 

j=0 

2n 

1 - [2n + 2](0)c + X][2« + 1 - ^](0 J> • • • J,i + l)c 

i=0 

2n 

= [n + l](0)c — ^['^ — ■j](0, 1, • • • , ^, ^ + l)c (easy arithmetic) 

i=0 

2n-l 

= [n + 1](0, • • • , 2n - l)s - n(0, • • • , 2n)s - ^ 



(Lemma 8) 



+n(0, • • • ,2n - l,2n,2n + 1, • • • ,4n)s (Lemma 9,r = 2n - 1) 
[n + 1](0, • • • , 2n - 1)b - n(0, • • • , 2n)B + nG2n 

-- [n + l](0,---,2n-l)B-n(0,---,2n)B (by (21)) 



(21) 



[n + 1] 



2n-2 



1 - 2n(0)B + ^ [2n - 1 - z](0, 1, • • • , z, z + 1)b 



-n 



;=o 

2n-l 



1 - [2n + 1](0)b + ^ [2n - i](0, 1, • • • , z, i + 1)^ 



i=0 



(dual of Lemma 
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2n-l 

= 1 — 71(0)5 + [n — 1 — i](0, 1, i, i + 1)b (easy arithmetic) 

2ri-2 

= 1 - n(0, ■ ■ ■ ,2^r^)A + [n- 1](0, • • • ,2n-l)A + 

i=l 

-n(0, • • • , 2n - 2, 2n - l,2n, • • • , 4n - 2) a (Lemma 10, r = 2n - 2) 

= c^-c^-i - nG2n-i (by (20)) 
which, upon adding nG2n on both sides, gives the claimed formula for (pCn ■ D 

As an example, let us compute the output distribution of C2. From A = \/^ CqX = \/^ X follows 
B = A = f\'\/ X. Using expansion calculus the reader may verify that 

G2n = G4 = (0, 1, 2, 3, 4, 5, 6, 7, 8)^ = pW + p'qf 

Similarly one gets 
Therefore 

ct>c, = </>Ci + 2(G4 - G3) 

= 3p3 + 3p4 _ 9p5 _^ 4^6 _^ 4^7 _ ]^Qp8 _^ 4^9 _^ gplO _ g^ll ^ 2pi2 

As to robustness, from the above representation of 4>C2 ^^"^ by using Corollary H] we obtain that 
C2 is lower robust of order 3 like U2- Similar to L2U2 discussed in 5.1, the upper robustness of 
C2 is not inherited from L2. Indeed, we have 

(j)(j^-l = q2(2pio _ 4p9 _ 2p^ + Ap'^ + ip'^ -p^ - 3p2 -2p-l) 

which implies that C2 is upper robust only of order 2. However, upper robustness is not con- 
stantly 2; these results were obtained from Theorem 12: 



n 


1 


2 


3 


4 


5 


6 


lower robustness 


2 


3 


4 


4 


5 


6 


upper robustness 


2 


2 


3 


3 


4 


4 



We mention that some handy closed formula for G2n and G2n-i is conjectured in [1, section 
4.5.6]. 
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6 Using some principle of exclusion as opposed to inclusion- 
exclusion 



As witnessed by 5.2, one can sometimes exploit symmetry to tame the inherent exponential 
complexity of inclusion-exclusion. However, without the possibility to clump together many 
identical terms, the number of summands in Lemma 7 is 2", which is infeasible already for 
n = 20 or so. 

In [14] on the other hand, some multi-purpose principle of exclusion (POE) is employed. When 
POE is aimed at calculating the output distribution of a stack filter S, a prerequisite is that 
the stack filteiB S be given as a disjunction of conjunctions Ki, i.e. in disjunctive normal form 
(DNF). The POE then starts off with the set Modi of all 0, 1-strings that satisfy Ki, then from 
Modi one excludes all 0, 1-strings that violate K2. This yields Mod2 ^ Modi, from which all 
0, 1-strings are excluded that violate K^, and so on. The feasibility of the POE hinges on the 
compact representation of the sets Modj. 

The details being given in [14], here we address the question of how one gets the DNF in the 
first place. Specifically we consider the frequent case that our stack filter 5 is a CDE (section 
4) whose structural elements are provided. Let us go in medias res by reworking S = \/^ /\^ 
of Example 1: 



{SX)o = Ao = ZoVZi (DNF) 

= (y_2Ay_iAyo) v(y_i AFo AFi) (blowup) 
= yoAF-i A(y_2vyi) (get cnf) 

= (Xo VXi VX2 VX3) A (X_i V Xo VXi VX2) (blowup) 

A ((X_2 V X_i V Xo V Xi) V (Xi V X2 V X3 V X4)) 
= (Xo V Xi V X2 V X3) A (X_i V Xo V Xi V X2) (condense) 

A (X_2 V X_i V Xo V Xi V X2 V X3 V X4) 
= (Xo V Xi V X2 V X3) A (X_i V Xo V Xi V X2) (condense further) 

= XoVXi VX2V(X3 AX_i) (get DNF) 

The last line is the sought DNF of S. It is obtained by starting with the DNF ZqW Zi. This gets 
"blown up" to a DNF in terms of 5^'s (using definition (13) of Zq and Zi). This DNF needs to be 
switcheqll to CNF (= conjunctive normal form). This in turn is blown up to a CNF in terms of 
Xj's. Usually the result can and must be condensed in obvious ways ("condense further" meant 
that only the inclusion-minimal index sets carry over). Like this one takes turns switching DNF's 



*More precisely, the positive Boolean function that underlies the stack filters must be given in DNF. 
^How one switches between DNF and CNF of a positive Boolean function (i.e. one without negated variables 
like here) is a well researched topic which we don't touch upon here. 
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with CNF's, and blowing up expressions. This is done as often as there are structural elements. 
As a bonus the so called rank selection probabilities rsp[i] are calculated. This is defined as the 
probability that the filter selects the i-th. smallest pixel in the ti;-element sliding window. For 
instance here w = 5 and rsp[l] = rsp[2] = rsp[3] = 0, rsj9[4] = 0.4, rsp[5] = 0.6. 

There exists a Mathematica Demonstration Project whereby the user provides the structural 
elements of any desired CDE S (also 2-dimensional) . From this the program first calculates the 
DNF, which then serves to calculate the output polynomial 4>s{p)- Alternatively the DNF of 
any stack filter (CDE or not) can be fed directly by the user. Albeit Wild's algorithm is multi- 
purpose, it managed to calculate 4>Cnip) up to n = 5 (and the result agreed with Butler's). 
Written out as CDE we have C5 = A^*^ ' ' ' V- "^^^ corresponding structural elements 
{0, 1, 2, 3, 4, 5}, {0, -1, • • • , -10}, {0, 1, • • • , 9} and so forth triggered the calculation of a DNF 
comprising a plentiful 12018 conjunctions (time: 168224 sec). From this (pCsip) was calculated 
in 45069 sec. Here it is: 

12x^ + 7x^ - 23x'^ + 19x^ - 130x^ + 194x^0 - 59x^^ - U2x^^ 
+460x^3 _ 787^14 ^ 7153,15 _ 7^,16 _ iQ^Q^n 

+1959x1^ - 2216x^^ + 208x^0 + STllx^i - 6748x^2 
+8412x^3 - 7587x2^ + 2023x2^ + 4680x^6 
-7903x27 + 8839x^8 - 13540x^9 + 30009x^0 
-51715x^1 + 50159x=^2 _ 7686x^3 - 51417x^4 
+78198x^5 _ 50589x^6 + 6900x^7 - 7680x3^ 
+56330x^9 - 86905x4° ^ 43710x^1 + 49540x^2 
-114680x^3 + 103390x^4 - 40555x^5 - 15370x^6 
+33955x47 - 25460x48 + 11790x49 - 3645x^0 
+740x5^ - 90x^''^ + 5x^3 



7 Conclusion 



The main result of this paper is the calculation of the output distribution of the LU LU -operators 
Cn and Fn- The basic technique for obtaining the output distribution in closed or recursive form 
may also apply to other highly symmetric erosion-dilation cascades; otherwise the algorithm of 
[14] does the job. Further we introduced the concept of robustness which is useful in character- 
izing the properties of a smoother in terms of its output distribution. 
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